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Abstract 

The Cellular Automaton (CA) modeling and simulation of solid dynamics is a long-standing 
difficult problem. In this paper we present a new two-dimensional CA model for solid dynamics. 
In this model the solid body is represented by a set of white and black particles alternatively 
positioned in the x- and y- directions. The force acting on each particle is represented by the 
linear summation of relative displacements of the nearest-neighboring particles. The key technique 
in this new model is the construction of eight coefficient matrices. Theoretical and numerical 
analyses show that the present model can be mathematically described by a conservative system. 
So, it works for elastic material. In the continuum limit the CA model recovers the well-known 
Navier equations. The coefficient matrices are related to the shear module and Poisson ratio of the 
material body. Compared with previous CA model for solid body, this model realizes the natural 
coupling of deformations in the x- and y- directions. Consequently, the wave phenomena related 
to the Poisson ratio effects are successfully recovered. This work advances significantly the CA 
modeling and simulation in the field of computational solid dynamics. 
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I. INTRODUCTION 



Cellular Automaton (CA) is a kind of discrete model for complex systems containing large 
numbers of simple components with local interactions [l-5]. It is charming that CA is able 
to present macroscopic complex phenomena with simple rules. The interesting fact hints 
that CA can provide a new conceptual framework, as well as an effective numerical tool, 
which can be used to explore the microscopic or mesoscopic mechanism of physical systems 
by simulating macroscopic behaviors. Further, we can capture more essential features of 
given physical systems by amending the rules. Therefore, the CA approach and the related 
modeling techniques are powerful methods to describe, simulate and understand physical 
systems. During the past decades CA models have been introduced and studied in various 
fields [lRlOj, such as the artificial life, traffic flow, quantum computation, snow avalanche, 
disease control and prevention, etc. 

n 

Historically, Frisch et al. [1| obtained the correct Navier-Stokes equations starting from 
the lattice gas CA and introduced CA approach into the field of fluid dynamics in 1980s, 
since then the CA approach has been particularly successful in this field. It is meaningful 
to mention that study of lattice gas CA has also inspired the appearance of the well-known 
Lattice Boltzmann (LB) method which has attracted extensive attention and has been suc- 
cessfully applied in various fields 

Besides the fields mentioned above, some researchers have developed the CA/LB approach 
to study viscoelastic fluid and solid. For example, Giraud, et al. proposed an LB model for 

n n 

viscoelastic fluid |13j. Mora, et al. presented their LB phononic lattice solid in Ref. |14J |. 
Fang, et al. proposed an LB model for photonic band gap materials [l5[]. Xiao studied the 
shock wave propagation in solids by LB method 16] . Chopard and Marconi et al. proposed 



a CA of solid which can be thought of as a grid of masses linked by springs [4], ll7|]. Based on 
the CA, they developed a kind of LB method for solids. Matic and Geltmacher developed a 
CA for modeling the mesoscale damage evolution |18| . 

However, in the field of solid mechanics, due to some unsolved problems, the CA approach 
has not been as successful as in the field of fluid dynamics. For instance, Chopard, et al. 
presented a very nice idea to recover the elastic behavior, but their CA still can not couple the 
motion in x direction with that in y direction. This limitation furthers to plague describing 
the Poisson ratio of solid materials and more practical behaviors. In order to simulate the 
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dynamics of elastic solid material, it is necessary to modify the solid CA of Chopard, et al. 
to couple the motions in different directions. This is the main aim of the present paper. 

The following parts of the paper are planned as follows. In the next section, we briefly 
review Chopard's CA, describe the improved model and perform numerical stability analysis 
for the new model. Some numerical examples are given in Sec. Ill to validate the model. 
Conclusions are summarized in Sec. IV. 

II. CELLULAR AUTOMATA MODEL FOR ELASTOMER 
A. Chopard's Cellular Automata 

A sketch of the two-dimensional particle set for Chopard's CA is shown in Fig. 1. In the 
case of translation motion, the automata rules are: (i) Divide the lattice into two sublattices 
(black/white) in a checkboard manner; (ii) Alternatively update all particle positions in 
each sublattice. Considering only the nearest neighboring inter-particle interactions gives 
the following rules: 




where b and w are the displacements of black particles and white particles relative to 
their initial equilibrium positions, respectively; Wj are the displacements of black particles's 
neighbor and vice versa; K is the number of a cell's neighbor; "+" denotes the next time 
step. In the case of rotation motion, the automata is quite difficult to describe the two- 
dimensional and three-dimensional behaviors. Through Chopard's CA above, we can see 
some mesoscopic mechanisms working in the framework of cellular automata modeling. In 
this section, we aim to obtain a CA that can represent the behaviors of elastomer ruled by 
Navier equation 

d 2 t M = (A + //) V(V • u) + (2) 
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FIG. 1: Grids for the two dimensional CA. 



Ax 
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FIG. 2: Directions of displacement vectors of a particle. 



B. Our new Cellular Automata for elastomer 



In the present work, we focus on the two dimension CA model, which is capable of coupling 
the motion in x-direction with that in y-direction and capable of representing the Poisson 
ratio of solid materials. We use the same particle set as shown in Fig. 1. The rules of CA are 
as follows: (i) Every particle updates its position according to the neighbors's current states; 
(ii) Subgrids (black and white) alternatively update particles' positions. Obviously, every 
particle has four neighbor particles in this model. The directions of every cell's neighbors 
are shown in Fig. 2, where Aj (i = 1,2,3,4) is the neighboring vector at the initial time. 
Now we can modify Eq. (1) to the following form 

b + = -b + B * w *> 

(3) 

w + = _ w + ^ Di b+, 



where Bj are the coefficient matrices of black particles' neighbors, Dj are the coefficient 
matrices of white particles' neighbors. It is noteworthy that Eq. (3) will degenerate into Eq. 
(1) when B; = D, = 1/2. 

From the basic evolution Eq. (JTJ), we can obtain the following volume ratio of consecutive 
two evolution steps in the phase space 

g({b+},{w+}) _ 
<9({b},{w}) " ' 

where {b} = {bj, i — 1,2, • • • } represents displacements of all black particles, {w} = 
{\Vj, i = 1,2,- ••} represents displacements of all white particles. It is clear that the CA 
system is conservative. 

Let Si = e AArV , Eq. (3) can be written as follows 



b + = -b + ^B^w, 

i 

w+ = - w + ^D^b + . 

Furthermore, if we let T = e XS9t , Eq. (4) can be rewritten as Eq. (5) 

(T+l)b=^B^w, 

i 

(T + l)w=^D^Tb. 

i 

To remove the quantity w by combining equations, we obtain 

(T + l) 2 b = BiSiDjSjTb. 

hi 



(4) 



(5) 



(6) 



To perform Taylor expansion with A on the both sides of Eq. (6) and compare the same 
order terms, we can obtain a series of equations. 
The zero-order term 

BD = 41, (7) 

where B = £\ Bj and D = ^D S . 
The lst-order term 



2(T + l)T'b = ^(B^D^T + BiSiDjS'jT + 'B i SiDjSjT')b. 



(8) 
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Eq. (8) can be rewritten as 

A5d t b = (BD + BD + BT>5d t )b, 

where 

i 
i 

Combining Eqs. (7) and (9), we have 

BD + BD = 0. 

Expanding Eq. (11), we obtain 

[(Bi - B 3 )D + B(Di + B 3 )]d x + (B 2 - B 4 )D + B(D 2 + D 4 )]d w = 0. 

The 2nd-order term 

2(T + l)T"b + 2T' 2 b = ^(BiS'jVjSjT + B^D^'T + B t S t B j S :j T")b 

+ 2 ^(BiS'iDjSfi + BiS'jDjSjT' + B^D^T')b. 

Eq. (13) can be rewritten as 

4{5d t ) 2 b + 2{5d t ) 2 b = [BD + BD + BB(5d t ) 2 ]b + 2(BD + BB5d t + BB5d t )b, 
where 

' B = ^B,(A l -V) 2 , 

i 

D = ^D J (A,-V) 2 . 

Eq. (14) can be changed into 

2(5d t ) 2 b = (BD + BD)b + 2BDb. 

Expanding Eq. (16), we have 

<9 2 b = ^R 2 [(B! + B 3 )D + B(D X + D 3 ) + 2(B 1 - B 3 )(D! - D 3 )] 
+ <9 2 [(B 2 + B 4 )D + B(D 2 + D 4 ) + 2(B 2 - B 4 )(D 2 - D 4 )] 
+ ^^[((Bx - B 3 )(D 2 - D 4 ) + (B 2 - B 4 )(Dx - D 3 )]}. 



In order to study the deformation of an elastic body, we select the Navier equation 



Ofh = c(V 2 h + c^VV • b 



4 + 4 



+ d] 



.2 | C l 



4 + 4 



+ d x d y 



r 2 

c 2 



(18) 



where C 2 = /z, c 2 = A + /i, A and ji are Lame coefficients. Physically, C\ is the sound speed 
of transverse wave, \J~4 + 4 i s the sound speed of longitudinal wave. 
Combining Eqs. (7), (12), (17) and (18), we get 

' D = B = 21, 
B 1 -B 3 + D 1 -D 3 = 0, 
B 2 - B 4 + D 2 - D 4 = 0, 



— [Bi + B 3 + Dx + D 3 + (Bi - B 3 )(Di - D 3 )] 



A 2 

— [B 2 + B 4 + D 2 + D 4 + (B 2 - B 4 )(D 2 - D 4 )] = 



— [(Bi - B 3 )(D 2 - D 4 ) + (B 2 - B 4 )(Di - D 3 ) 



4 + 4 



(19) 



4 + 4 



r 2 



For the sake of convenience, we let ^ = 1, thus 4 an d 4 are dimensionless. Consequently, 
Eq. (19) is simplified to 

D = B = 21, 
Bi-Bg + D 1 -D 3 = 0, 
B 2 - B 4 + D 2 - D 4 = 0, 

Bi + B 3 + Dx + D 3 + (Bi - B 3 )(Di - D 3 ) = 
B 2 + B 4 + D 2 + D 4 + (B 2 - B 4 )(D 2 - D 4 ) = 
(Bx - B 3 )(D 2 - D 4 ) + (B 2 - B 4 )(D! - D 3 ) 



4 + 4 



(20) 



4 + 4 



r 2 

c 2 
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Now, we obtain the discrete model described by Eq. (20). By selecting a series of appropriate 
matrices Bi and A, we will obtain a CA model for elastic body. 

To seek such matrices, we introduce another two matrices P and Q that are defined as 
in Eq. (21) 

' Bi - B 3 = D 3 - D x = P, 
B 2 - B 4 = D 4 - D 2 = Q. 

Thus, Eq. (20) can be rewritten as 



, 2cj + c 2 2 




2c\ + 4 



and Eq. (22) can be changed into 
41 - (P + Q) 2 = 



2c 2 + 4 



2c{ + 4 





< 4\ 


H 













Here we select the forms of P and Q as follows 

P = p l + Pi • <7i, 

^ Q = q l + q 1 ■ (Ti, 

where 

o"i = 




is the first Pauli matrix. So we have 

P 2 + Q 2 = (Po + Qo +Pl + & + 2(poPi + gb9i)<Ti, 
PQ + QP = (2p g + 2pigi)I + 2{p oqi + q Pi)(Ji. 

Therefore, Eq. (23) can be changed into algebraic equations 

pl + ql+pl + ql = 4-2c 2 -c 2 , 

PoPi + <?o<?i = 0, 

poqo+piqi = o, 
- 2(p gi + qopi) = 4- 



(21) 



(22) 



(23) 



(24) 



(25) 



(26) 



(27) 



After choosing appropriate c\ and c 2 , we can obtain the values of po, go, Pi, an d q±, namely 
we can determine P and Q by solving Eq. (27). 

Now we can determine the series of matrices Bj and D» by solving the following equation 

B x + B 2 + B 3 + B 4 = 21, 
D x + D 2 + D 3 + D 4 = 21, 
B x - B 3 = P, 

Dl - D 3 = -P, (28) 
B 2 — B 4 = Q, 
D 2 -D 4 = -Q, 
B 1 + B 3 + D 1 + D 3 = M + P 2 , 



where 



M = 



cj + 4 



(29) 



Need to emphasize again that c\ and c 2 have the ability to describe the elementary prop- 
erties of elastomer. Therefore, the CA model is able to describe the elementary properties 
of an elastic body. 

Through selecting a series of po, p±, qo, and q\ as follows 



Po = -qi = -cj, 



pi = go = yi- \c\ - \c\. 

we can obtain a CA model with the following parameter matrices 



B x 


= D 3 = 


i(M + 2P + P 2 ), 


B 2 


= D 4 = 


i(4I-M-P 2 + 2Q), 


B 3 


= D 1 = 


i(M-2P + P 2 ), 


B 4 


= D 2 = 


i(4I-M-P 2 -2Q). 



(30) 



(31) 



Especially, when c\ = y/2 and c 2 = 0, the matrices Bj = Dj = |I. Thus, the model degrades 
into Chopard's model. 
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C. Stability analysis of the Cellular Automata model 



iW 



The numerical scheme of this CA model can be written as follows 

b+ = -b + ^B^w, 

i 

W+ = -W + D A (-b + B i S i 

j i 

We change Eq. (32) into matrix form 

-I E,b^ 

Performing Fourier transform to Eq. (33), we have 



(32) 









l w) 



(33) 





\=G k \ 















(34) 



where 



Gk — 



—I Ei BjSfci 

E, Dj-S'fcj Ei, 7 - Dr^/ lirS.-; - I 



(35) 



Here BjSfcj = Bj-e lk J ', T^jSkj = Dj-e* k is the wave vector, and Aj is the displacement 
direction vector. 

Selecting appropriate matrices Bj and Dj, we solve the eigenvalues of 



= 0, 



(36) 



-I - AI EjBjS^ 
Ej Dj-S'fcj Ej,j D i S' fcj BjS' fc j - I - AI 

and let |Aj| = 1, thus we can obtain a conservative linear CA model. The model described 
by Eq. (3), Eq. (24) and Eqs. (29)-(31) is a conservative linear system. Namely, this model 
presents the ideal wave process in linear elastomer. 



III. NUMERICAL EXAMPLES 

In this section, several numerical examples are used to validate the CA model. The 
former two describe the translational motion in elastomer, and the latter one simulates the 
single-mode motion in elastic body. 
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FIG. 3: CA simulations of the translational motion within the region. 
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FIG. 4: CA simulations of the translational motion with bounce to boundaries. 
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A. Examples of translational motion 

In the former two examples, we select a square area 20 x 20 as the computing domain. 
The size of mass is 6 x 6. Here we set C\ = 0.4, c 2 = 0.58. Boundary conditions (BCs) are 
implemented as follows: (i) For the solid to bounce, any x or y motion over the boundary 
is simply forbidden; (ii) Free BCs are imposed on object which has some ghost particles on 
the boundary. 

When the distance in x or y direction between the neighboring two particles equals to a 
constant A , the mass is static. Here A = 1.5. The initial disturbance is in the following 
form: x'(3, 2) = x(3, 2) — 0.2A , as shown in Fig. 3. It is noteworthy to mention that the 
initial disturbance is only added in x direction. 

There are two numerical examples about translational motion which share the same initial 
state. Example (I) presents the evolution of the CA model within the region, as shown in 
Fig. 3. Example (II) presents the translational motion with bounce to boundaries, as shown 
in Fig. 4. In Example (I), although the disturbance is only added in x direction, from Fig. 
3 we can see that particles of the mass present motions in y direction. The couple between 
different directions is illustrated more clearly in Example (II). As a result of the contact with 
the boundary, not only the particles have motions in x direction, but the whole mass has 
motion in y direction (see Fig. 4). Therefore, we conclude that the CA model can couple 
the motion in x direction with the motion in y direction. 

B. Example of single- mode motion 

To further validate the model, an example of single-mode motion is illustrated here. In 
this example we select a square area 20 x 20 as the computing domain. The size of object 
is 8 x 8. In this case, we set c\ = 0.5 and c 2 = 1.0. 

In this example, Ao = 1.0 when the mass is static. To obtain the single- mode motion, 
we select the eigenvalue A and the corresponding eigenvectors h k0 and w fc0 as follows 

A = 0.75 + ^0.661438, (37) 
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FIG. 5: Comparisons between CA simulations and the exact solutions. 



14 



\ a-22 J 



( 0.176777 + ^0.467707 \ 
-0.176777-^0.467707 
-0.5 
0.5 



(38) 



The initial displacements of the particles in mass are given by the following equation 



^ box \ 




/ 


( an ^ 


\ 




= Re 






ikr 








021 




\W y J 




V 


V °22 / 


/ 



(39) 



where k are the wave vectors, and r are the situation vectors of every particle before the 
mass is disturbed. From the view of analytic solution, the motion of the mass will be ruled 
by the following equation 

/ b x (t) \ 1 1 a n \ \ 



b y {t) 
w x (t) 

\Wy(t)J 



Re 



V 



Ol2 
021 

\ a 22 ) 



i(k-r+tot) 



(40) 



/ 



where oj = — i In A. 

The numerical scheme is ruled by Eq. (3), Eq. (24) and Eqs. (29)-(31). The initial 
configuration is given by Eq. (39), as shown in Fig. 5. The BCs are simpler than the above 
two examples. Here we need not to consider the contact between mass and boundary of the 
computing domain. BCs of the mass are periodic. Comparisons between CA solutions and 
the exact solutions are illustrated in Fig. 5. It is clear that the two sets of results have a 
satisfying agreement. Further, we draw the conclusion that the CA model has the ability to 
present qualitatively wave propagation in elastomer. 

It is necessary to point out that although the proposed model recovers the Possion ratio 
effects, compared with the deformation uniformity and periodicity of the true elastomer, 
the CA simulation results (see Figs. 3-5) still exist some gaps. These gaps firstly arise 
from the lack of introducing appropriate rules describing rotaion and elastic-plastic of solids 
in the CA model. Secondly, it should be emphasized that what the CA model concerns 
are different with other macroscopic numerical methods (e.g., finite-element method). The 
CA approach focuses on understanding of the physical mechanism behind the macroscopic 
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physical phenomena, but not the comparisons with specific experimental data. Indeed, the 
CA approach can be regarded as coarse and large time step sampling results of the elastic 
body during its movement. Therefore, there are some differences between the CA simulation 
results and the realistic elastomer. 

Meanwhile, we should note that the present work has already pointed out the way how to 
introduce the above-mentioned physical characteristics into the CA model. In other words, 
the new CA model can be further modified to describe more realistic dynamical behaviors 
in solid materials. For example, via relating the coefficient matrices to the displacement 
field, the CA model can be used to describe the nonlinear elastic phenomena; via changing 
the coefficients of b and w from the fixed value, —1, the CA model can be used to describe 
the inner dissipation in the plastic behavior of the solid body; via breaking the neighboring 
relations of particles, the CA model can be used to describe damage and fracture phenomena. 
Generally speaking, since the mesoscopic collision rules of the CA model and the macroscopic 
behaviors of solid exhibit complex nonlinear relationships, looking for mesoscopic rules which 
can be used to describe macroscopic behaviors (rotation, elastic-plastic, etc.) of solid is still 
an open problem. Recently, the progress has been rather slow and there is rare report on 
this issue. So, from this point of view, this paper provides us with an effective way to 
propose more powerful CA models and study dynamics of solids via this more essential and 
fundamental approach. This is another contribution of the present work. 



IV. CONCLUSIONS 



A new two-dimensional cellular automaton model for solid dynamics is proposed. The 
solid body is represented by a set of white and black particles alternatively positioned in the 
x- and y- directions. The force acting on each particle is represented by the linear summation 
of relative displacements of the nearest-neighboring particles. The key technique in this new 
model is the construction of eight coefficient matrices. Theoretical and numerical analyses 
show that the present model can be mathematically described by a conservative system. 
So, it is a model for elastic material. In the continuum limit the CA model recovers the 
well-known Navier equations. The components of the coefficient matrices can be related to 
the elastic parameters of material, such as the Lame coefficients or the shear module and 
Poisson ratio. 
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Chopard et al. introduced the CA model into the field of solid dynamics for the first 
time. Even though in the form of two dimensions, the original CA model by Chopard et al. 
is in fact quasi-one dimensional and can only describe some simple wave phenomena. The 
new model presented in this paper realizes the natural coupling of deformations in the x- 
and y- directions. Consequently, the wave phenomena related to the Poisson ratio effects are 
successfully recovered. When the eight coefficient matrices are all diagonal, the new model 
goes back to the original CA model by Chopard et al. This work advances significantly the 
CA modeling and simulation in the field of computational solid dynamics. 
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